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Abstract 

A power series expansion in the damping parameter € of the limit cycle 
U(t;€ ) of the free van der Pol equation U+e[U-lJU+U = 0 is 

constructed and analyzed. Coefficients in the expansion are computed up to 
0(€^) in exact rational arithmetic using the symbolic manipulation system 
MACSYMA and up to 0(e^^) using a FORTRAN program. The series is analyzed 
using Fade' approximants. The convergence of the series for the maximum 

amplitude of the limit cycle is limited by two pair of complex conjugate 
singularities in the complex G-plane. A new expansion parameter is introduced 

which maps these singularities to infinity and leads to a new expansion for 

the amplitude which converges for all real values of €. Amplitudes computed 
from this transformed series agree very well with reported numerical and 

asymptotic results. For the limit cycle itself, convergence of the series 

expansion is limited by three pair of complex conjugate branch point 
singularities. Two pair remain fixed throughout the cycle, and correspond to 
the singularities found in the maximum amplitude series, while the third pair 
moves in the € -plane as a function of t from one of the fixed pairs to the 
other. The limit cycle series is transformed using a new expansion parameter, 
which leads to a new series that converges for larger values of € . 
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1) Introduction 

We wish to study the limit cycle U(t;e) of the free van der Pol equation 
(1.1) U + e(U 2 -l)U + U = 0, 

where the dots represent differentiation with respect to the independent 
variable t . We shall do this by constructing and analyzing the power series 
expansion of U in the damping parameter e . We first wish to determine some 
of the analytical structure of U as a function of e and, in particular, to 
determine the locations of the singularities of U in the complex e-plane 
which are closest to the origin and, hence, limit the radius of convergence 
of the power series solution. Once the locations of these singularities 
have been determined, a new expansion of U can be constructed which will 
converge for larger values of e. 

This work is a companion paper to some work reported previously by 
Andersen and Geer [1] (henceforth referred to as I) on the power series 
expansions of the frequency and period of the limit cycle of the van der 
Pol equation. In I, the power series expansions of the frequency v of the 
limit cycle U, as well as U itself, were computed up to a high number of 
terms. However, in I, only the series expansion of v was analyzed in detail. 
In particular, it was found that the convergence of the series for v was 
limited by a pair of complex conjugate branch point singularities in the 
complex e 2 -plane. These branch points, which were located at e 2 = Re - '*’^ 
with R = 3.42 and 3 =1.7925, appear to be the only singularities of v in 
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the finite part of the e -plane. Hence, when a new expansion parameter was in- 
troduced which mapped these singularities to infinity, the resulting expansion 
converged for all real values of e. The values of the period computed from a 
completed form of this new expansion compared very well with reported numeri- 
cal results, as well as with the asymptotic formula for the period valid for 
large values of e . 

In section 2 below, the problem of determining the limit cycle is formu- 
lated and the method we used to compute the power series expansions of V and 
U is described briefly. In section 3, the series for the maximum amplitude of 
U, corresponding to U at t=0, is analyzed using Pade" approximants . It is found 
that the convergence of the series is limited by the presence of two pair of 

complex conjugate branch point singularities in the complex e-plane, which 

2 

correspond to the square roots of the singularities in the e -plane found for 
the frequency series. One pair of these singularities lies in the first and 
fourth quadrant, while the other pair lies in the second and third quadrant. 

A new expansion parameter is introduced in section 4 which leads to an ex- 
pansion for the amplitude that converges for all real values of e. Values of 
the amplitude computed from this series, and a modified version of it, are 
computed and found to compare very well with reported numerical results and 
also with values computed from the asymptotic formula for the amplitude valid 
for large e . 

The limit cycle itself is analyzed in section 5 using Pade approximants 
for values of t between 0 and T, where T is the period of the limit cycle. 

We now find that there are three pair of complex conjugate branch point singu- 
larities which lie about equidistant from the origin. For all values of t, 
two pair of these singularities remain fixed in the e-plane at the same 
locations as the singularities for t=0 (i.e., for the amplitude series). At 
t=0, the location of the third pair of singularities coincides with the loca- 
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tipn of the fixed singularities in the second and third quadrants. However, as 
t increases from zero, these singularities begin to move away from the fixed 
locations in the second and third quadrants toward the fixed locations in the 

first and fourth quadrants. In particular, they cross the imaginary axis in 

the complex e-plane when t=T/4 and arrive at the fixed locations in the first 
and fourth quadrants at t=T/2. This motion of the singularities repeats itself 
for T/2 < t < T. Knowledge of the location of these singularities suggests a 
new expansion parameter and a transformation of the original series which will 
converge for larger values of e. This is done in section 6, while our results 
are discussed in section 7. 

2) Power Series Expansion in e 

We construct the power series expansion of the limit cycle of the free van 
der Pol equation (1.1) by first making the change of variables 

(2.1) x = vt, 

where v =v(e) = 2ir/T(e) is the (unknown) frequency of the limit cycle and T(e) 
is its (unknown) period. We then let 

(2.2) u(x;e) = U (t;e) , 

where u is now periodic in x with period 2 tt. In terms of u, equation (1.1) becomes 

(2.3) V 2 u +ev(u 2 -l)u + u = 0, 

where the dots now denote differentiation with respect to x. In addition, we 
impose the phase condition that u have a maximum at x=0, i.e., 

(2.4) u(0,e) = 0, u(0,e) = A(e) > 0. 

We now look for solutions for v(e),^ u(x;e), and A(e) in the form 

(2.5) V (e) = 1 + Z v .e J , 

1 = 1 J 

(2.6) u(x;e) = £ u.(x)e , 

j=0 

(2.7) A(e) = ? a.e j . 

j=0 3 
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) 


Here, each Vj and is an unknown constant, while the unknown functions u^ (x) 
must be periodic with period 2tt. Substituting (2.5) and (2.6) into (2.3) and 
then (2.6) arid (2.7) into (2.4), we are led to a sequence of linear problems 
from which the u^ can be determined recursively. In particular, demanding 
that each of the u^ is periodic in x leads to the unique determination of the 
constants V ^ and a^ (see I, section 2, for details). In this way we find 


v 


2i+l 


® * 1 


while the u^ (x) have the form 


2i 

u 9 . = T. u 9 . . cos ( (2k+l)x) , 
k=0 5 


2i+l 

u 2i+l = . u 2i+l,k 
k=0 


sin((2k+l)x) , i < 0,1,2,..., 


where the u. , are certain constants. 

J ,k 

Through the use of the symbolic manipulation system MACSYMA we have de- 
termined the coefficients in the expansions of both v(e) and u(x,e) up to 
24 

0 (e ) in exact rational arithmetic. The exact rational expressions for the 
coefficients v ^ are reported in I, while the first few terms in the expansions 
of A(e) and u are 


. .,12 1033 4 . 

(2.8) A - 2 + — e - 55296 o€ + 


1019689 


55738368000 

9835512276689 8 

157315969843200000^ 

58533181813182 818069 
7 32614 17 89 209 886 7 20000000 6 


10 


12 

+ 0 ( 6 ^) , 


(2.9) u = 2cosx + (ysinx - 7^sin3x)e + (--jrcosx + A cos3x - cos5x)e^ 

4 s 8 lb yb 


, / 7 . . 21 . „ 

+ (-Trrrsmx + -^rsm3x 


256 


256" 


-^^sin5x + ^ygSin7x)e^ 


(73 

+ 12288 COSX 


47 „ , 1085 

rCOS3x + ■ —> 7 O cos5x 


1536 


27648 


2149 

110592 


cos7x + 


^ cos9x)e^ 


20480 


+ 0(c J ). 



- 6 - 


We then constructed a FORTRAN program (using floating-point arithmetic) which 

X6 3 

evaluated the various coefficients up to terms which are 0(e J ). (See I, 

section 2, for some details which expedited the computations.) In particular 
the coefficients a^ ,a ^ , . . . » a ^g computed from this program are listed in floating 
point form in Table I. 

Using some of the same ideas as in I for the frequency series, we now wish 
to investigate the convergence of the series (2.6) and (2.7). In particular, 
we wish to determine the location of the singularities of u(x,e) in the com- 
plex e-plane which are closest to the origin, and, hence, which limit the 
convergence of these series. 


3) Analysis of the Amplitude Series 

\ 

We begin our analysis of the limit cycle by considering the series for the 

maximum amplitude A(e). Using the phase condition (2.4), it follows that 

A(e) = u(0;e) , or 

00 2i 

(3.1) A(e)=E a^.e , where a^ . = v^.(0). 

4 —A 3 J J 


Values of A(c) determined by (3.1) for e between 0 and about 2 are indicated by 

the dotted line in Fig. 1. Clearly, the series diverges for e close to 2. The 

root test on the coefficients a„ . provides the estimate that the radius of con- 

^3 

vergence R q of the series (3.1) is approximately 1.86. 

In order to investigate the convergence of the amplitude series more 
closely, a sequence of [N/N] Pade approximants were constructed from the Taylor 
series expansion (3.1) for N = 20 ,24 ,28, . . . ,64 . The location of the zeros and 
poles of the approximants in the complex e-plane were determined for each value 
of N. Since A(e) is a real function of e, any of the zeros or poles which did 
not lie on the real axis had to appear in complex conjugate pairs. In Fig. 2, 
we have shown the location of the zeros and poles for the [48/48] Pade 
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npproximant which lie in the upper Im I f (-plane. The pattern of nearly over- 
lapping zeros and poles in Fig. 2 is similar to the pattern observed in the 
analysis of the series for v(<) in I and indicates the presence of two pair 
of branch point singularities located at 

(3.2) c = ±R e ±li3 °, R * 1.85, 3 £ 0.8770. 

o o o 

Thus, we see that the convergence of the series (3.1) is limited by the 

presence of two pair of complex conjugate singularities located (approximately) 

at the points indicated by (3.2). These singularities also provide an estimate 

of the radius of convergence R , which is consistent with the estimate obtained 

o 

from the root test. We note that these singularities are very close to (most 

probably equal to) the singularities found in the series expansion of v(e) 

2 ~ +±Q 

in I. (In I, the singularities of v(e) were found in the form e = Re~ , or 
e = ±R^e ±±e/2 , with = (3.42)^ = 1.849 and 8/2 = 1.7925/2 = 0.89625. These 
values are close to those given in (3.2).) In Fig. 1, we have indicated by 
the solid line the amplitude predicted by the [48/48] Pade approximant for 
values of e between 0 and 9. • We note that the Pade approximant agrees well 
with the series (3.1), where it converges, and also agrees well with some 
values for the amplitude computed by purely numerical methods. (These 
numerical methods will be discussed in more detail in the next section) . 


4 ) Transformation of the Amplitude Series 

We can now use the information obtained in the previous section to in- 
troduce a new expansion parameter which will lead to a new series representa- 
tion for A(e), which will converge for larger values of e. To do this, we 
introduce a new expansion parameter 6 = 6(e), defined by the transformation 
(4.1) 6(e) = 


{R 4 -2R 2 e 2 cos28 


o o 

where R q and 8 Q are given (approximately) by (3.2). The transformation (4.1) 
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has the properties that the origin remains fixed (i.e., 6(0) = 0), the singu- 

+ 4 g 

larities IR^e - 0 in the e-plane are mapped to infinity in the complex 6-plane 
(i.e., |6(±R o e“ p °)[ = 00 ) , and the real axis in the e-plane is mapped onto the 
portion of the real axis in the 6-plane between -1 and +1 (in particular, as 
€•> + «>, 6*> + l). In addition, it is a simple task to invert the transfor- 
mation (4.1) to express e as a function of 6. 

Using the definition (4.1), we recast the series (3.1) for A(e) into a 
series in powers of 6 of the form 

(4.2) A(e) = S a e 2j = S a 5 2:j 
j=0 j-0 3 


where the new coefficients a can be expressed in terms of the coefficients 

^ J 

a„ , and the transf ormation (4.1) in a straightforward manner. Since the 
^ J 

±ig 

singularities at ±R o e ° appear to be the only singularities of A(e) in the 
finite part of the e-plane, we expect that the series (4.2) will converge for 
all values of 6 with 1 6 1 <1 . On the other hand, since A(e) has an essential 
singularity at e equal to infinity (see [2], for example), we would not ex- 
pect the series (4.2) to converge for values of 6 with | 6 | >1 . 

To investigate this point quantitatively, we performed a Pade analysis 
on the series (4.2) using the coefficients . The zeros and poles of the 
Pade approximants again formed patterns indicating branch point singularities 
(in the complex e-plane), this time lying on the real and imaginary axes at 
distances from the origin slightly greater than one. (For example, for the 
[48/48] approximant, the closest zero-pole pair was on the imaginary axis 
at a distance of 1.015 from the origin.) Thus, it appears that the radius 
of convergence for the series A(6) is nearly one. The root test on the co- 
efficients gives a result consistent with the value of one. We have 
used the series (4.2) to compute values of A for selected values of e between 
1 and 100, which correspond to values of e for which good numerical results 
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are available (see, e.g., Zonneveld [7]). The results are presented in the 
third and fourth columns of Table IT. and show that the series (4.2) agrees 
well with the numerical results, with a maximum error of about 0.79% at 
e = 100. In Figure 3, we have plotted A(6) using (4.2) (dotted line) for 

5 close to one. Note that the vertical scale in Fig. 3 is 100 (A-2), and, 
hence, the agreement between A(6) and the numerical results is better than 
might be indicated at first glance. 

We can now use some information about the asymptotic behavior of the 
amplitude as £ °° to improve the convergence of our series for A(6) near 

6 = 1 (i.e., for large values of e). As £ -*■ 00 , we know (see [5], for ex- 
ample) that 

(4.3) A = 2 - |e~ 4/3 - f<3e-l + 21og2 - 81og3)e“ 2 + 0(e _8/3 ), 

where a =-2.3381 is the highest zero of the Airy function and g = 0.1723. 
Since 6 = 1 + Ce + 0(e 4 ), as £-><», the form of (4.3) suggests that we 
write A(6) as 


(4.4) A(6) = 2+ (l-6 2 ) 2/3 Z b 6 2j , 

j=0 3 

where the constants b„. can be determined from the coefficients iL.. We note 

2j 2j 

that the representation (4.4) has the properties that A = 2 when 6 = 0 or 
6 = 1 (i.e. , when e = 0 or e = °°) and also that, as 6 -* 1, or e -> °°, the 
asymptotic form of (4.4) agrees with the form of the first two terms in (4.3). 

Amplitudes computed from (4.4) are listed in the fifth column of Table 
II for values of e between 1 and 100. We see that, for values of e greater 
than 8, the series (4.4) agrees more closely with the numerical results than 
does the series (4.2). In Table III \<re have compared values of the amplitude 
computed from the series (4.2) and (4.4) with the four term asymptotic formu- 
la (4.3) for A(e). The differences between the amplitudes computed from (4.2) 
and the asymptotic formula (4.3) steadily increase as £ increases, approaching 
a maximum error of about 0.85% at infinity. However, when (4.4) is used, 
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these differences steadily decrease, the agreement being perfect (by design) 
at 6 = 1 (i.e., infinite e). In Figure 3, we have also plotted the values of 
100(A-2) using (4.4). Again, the close agreement as 6 -> 1 among the values 
for the amplitude computed from (4.4), (4.3), and the numerical results can 
be clearly seen. 


5) Analysis of the Limit Cycle Series 

We now investigate the convergence of the limit cycle series (2.6) for 
values of x between 0 and it. (Since each u^(x-Hr) = -u^(x) for all x, the be- 
havior of the convergence of the series for x between it and 2tt will be the 
same as the behavior for x between 0 and ir) . In Figure 4, we show the phase 
plane plots of the limit cycle, using the first 100 terms in (2.6) for e = 0, 
1.0, 1.5, and 1.6. The plots for e = 1.0 and 1.5 are smooth and agree very 
well with similar plots computed by purely numerical means. However, for 
e = 1.6, we notice that, while most of the cycle is smooth and, in fact, agrees 
very well with the true cycle, portions of the curve in the second and fourth 
quadrant have developed oscillations which are not present in the actual limit 
cycle. These oscillations become more pronounced as e is slowly increased 
above 1.6, although certain portions of the curve still remain smooth. For 
e = 1.8, however, the limit cycle computed from (2.6) is completely meaningless. 

The non-uniform behavior of the convergence of the series representation 
of the limit cycle suggests that its radius of convergence may vary with x. 

In particular, the portions of the phase plane plots where the (artificial) 
oscillations appear correspond to values of x near tt/2 and 3tt/2, This suggests 
that the radius of convergence of (2.6) may be somewhat smaller for values of 
x near ir/2 and 3tt/ 2 than for other values of x. 

In order to investigate this behavior more carefully, we performed a Pade 


analysis on the coefficients u^ (x) for values of x between 0 and tt, in increments 
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of it / 50 . For values of x between 0 and about tt/4, the analysis indicated the 
presence of only the same two pair of complex conjugate branch points which 
were present for the amplitude series, i.e., for x = 0 (see Figure 2). In 
Figure 5, (a) and (b) , we show the location of the zeros and poles of the 
[24/24] Pade approximant which lie in the upper half complex plane for x = tt/50 
and x = 13tt/50. As x was increased to values greater than tt/4, however, 
something quite different and very interesting happened. It appears that, as 
x increases above tt/4, a "new" singularity leaves its position in the second 
quadrant and moves toward the imaginary axis. In Figures 5 (c) and (d) we 
show the zeros and poles of the [24/24] Pade approximant for x = 18tt/50 and 
x = 24tt/50, which clearly indicate the presence of this "moving singularity", 
while also indicating that there is still a singularity which remains fixed 
in the second quadrant. At x = tt/2, the moving singularity crosses the 
imaginary axis. As x increases above tt/2, the moving singularity enters the 
first quadrant, and moves toward the position of the "fixed singularity" in 
that quadrant. At about x = 3ir/4, it appears to coalesce (at least approxi- 
mately) with the fixed singularity and remains there for x between 3tt/4 and 
TT. As x varies between tt and 2tt, the phenomena described above exactly re- 
peats itself, with a singularity again moving from the second to the first 
quadrant. Of course, a similar phenomenon is happening simultaneously in the 
lower half plane. 

Thus, it appears that the limit cycle series (2.6) has three pairs of com- 
plex conjugate singularities in the complex e-plane, two pairs remaining fixed 
as x varies, while the third pair moves as a function of x. We have illustrated 
this schematically in Figure 6, indicating in the upper half complex e-plane the 
location of the fixed singularities and also the approximate location of the 
moving singularity for various values of x. 

From a knowledge of the location of these singularities, we can estimate 

\ 
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the radius of convergence R of the series (2.6). In particular, tjhe radius of 

) 

convergence will be just the distance from the origin to the closest singulari- 

3tt 

ty. For 0 < x < ir/4 and 4 < x < tt, when the moving singularity coincides (at 
least approximately) with the fixed singularities, the Pade analysis gives 
R = 1.85. However, for tt/4 < x < 3tt/4, the moving singularity becomes the 
dominant singularity, moving closer to the origin than the fixed singularities. 
Thus, R decreases for tt/4 < x < tt/2, reaching a minimum of about 1.65 at 
x = tt/2, and then increases to about 1.85 at x = 3tt/4. (We compared the values 
for R we obtained from the Pade analysis with values of R computed using the 
root test and found the results to be consistent.) We have plotted R as a 
function of x in Figure 7 (a). As mentioned earlier, the smaller radius of 
convergence of the series (2.6) for x near tt/ 2 accounts for the presence of 
the oscillations in certain portions of the phase plane plots in Figure 4. 

In particular, for e = 1.6 (Figure 4 (d)), the oscillations occur for values 
of x near tt/2 (and 3n/2), indicating that this value of e is close to the 
radius of convergence in this region. However, for x not near tt/2 or 3tt/2, 
g = 1.6 is far enough from the radius of convergence to render the first 100 
terms in (2.6) (used to calculate the plots in Figure 4) a good approximation 
to the true solution. 

In Figure 7 (b) we have plotted the angle 8 as a function of x. Here 
8 is the angle which the moving singularity makes with the positive real axis. 
Although we have not been able to determine with any great certainty what the 
moving singularity is "doing" for 0 ^ x < tt/4, there is some evidence that it 
may possibly be spiraling out from its location at x = 0 . One piece of evidence 
to support this conjecture is the little "bump" in Figure 7 (b) near x = tt/4, 
as well as smaller oscillations for x < tt/4 which are too small to show up on 
the figure. Similar, but less pronounced, oscillations occur near x = tt/4 in 
the plot of R(x) . Of course, similar oscillations occur near x = 3tt/4, 
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indicating that the singularity may spiral into its location at x = tt in the 
first quadrant. However, since we are dealing with very small oscillations 
and since we have two pairs of singularities coalescing, our numerical results 
are not accurate or stable enough to make any definite statement about this 
conjecture. 

6 ) Transformation of the Limit Cycle Series 

Using our knowledge of the location of singularities of the limit cycle 

in the complex e-plane, it is possible to introduce a new expansion parameter 

which will allow us to recast the original series (2.6) into a form which will 

converge for larger values of e. To do this, we let e = ±R Q e~ ° be the loca- 

±ig 

tions of the fixed singularities and e = Re be the locations of the moving 
singularities. Here R q = 1.85 and 8 q = 0.8970, while R and 8 are functions of 
x, as indicated in Figure 7. If we now define 

(6.1) w = w(e) e , 

{(R 4 - 2R 2 cos 28 e 2 + e 4 ) (R 2 - 2Rcos 8 + e 2 )} 1/6 
o o o 

we see that w(0) =0, |w| -> 1 as |e| -* °°, and |w| ->-<»as e approaches any one 
of the six singularities of u. Using (6.1) we can recast (2.6) into the form 

CO ,00 

(6.2) u(x;e) = E u.(x)e 3 = E u.(x)w 3 , 

j=0 3 j=0 3 

where the functions u^. can be determined in a straightforward manner. We 
performed a root test on the coefficients u_. for values of x between 0 and 
Tr and found in each case an estimated radius of convergence very close to 
one. We also performed a Pade analysis of the coefficients u and found 
again a constant radius of convergence very close to one. In particular, 
for each value of x, the Pade analysis indicated singularities of u in the 
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complex w-plane very near w = ±1. These singularities are related, of course, 
to the essential singularity of u at e equal to infinity. 

To obtain some phase plane plots using the transformed series (6.2), we 
found it more convenient to use a transformation which is simpler (but less 
accurate) than (6.1). (For practical computations, the use of (6.1) requires 
the recomputation of the series coefficients u_. for each value of x used, 
since R and 3, and hence w(e), depend on x.) Instead of (6.1), we used the 
transformation 


(6.3) 


w = w(e) = 


{(R 4 - 2R 2 cos 23 e 2 + e 4 )(R 2 + e 2 )} 1/6 
o o o m 


where R m = min R(x) = 1.65. Here w has the same mapping properties as w, except 
that it maps the moving singularity to infinity only when it is closest to the 
origin, i.e., when x = tt/2. Thus, the region in the e-plane containing the 
singularities is mapped out away from the origin, but not all the way to infini- 
ty. The advantage of (6.3), of course, is that it is independent of x and needs 
to be inverted only once during the computations to express e as a function 
of w. 

Using (6.3), we recast the series (2.6) into a form similar to (6.2), with 

w replaced by w. A Pade analysis of this transformed series revealed that it 

/ 

has a minimum radius of convergence of about 0.953, which corresponds to a 
radius of convergence in the e-plane of about 3.97. In Figure 8, we show the 
phase plane plots of u using this transformed series for e = 1.65 and e = 3.5. 


7) Discussion 

The method of analysis we have presented here is an interesting combination 
of the methods of regular perturbation analysis and classical analysis, with 
both symbolic and numerical computation used to carry out the details of the 
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methods. As we mentioned in I, the importance of the role of symbolic computa- 
tions in this analysis should not be underestimated. Among other advantages, it 
allowed us to find certain patterns and forms of the solutions which greatly 
facilitated the construction of an efficient numerical FORTRAN program to carry 
the computations further. Also, the coefficients in the various expansions 
obtained by symbolic computation, being in exact rational number form, allowed 
us not only to check the results of our FORTRAN program but also to estimate 
the magnitude and effects of roundoff errors in our numerical computations. 

As far as the practical application of our results is concerned, it is 
interesting to note that we have been able to present formulas for the fre- 
quency and limit cycle function, valid for large values of e, by knowing only 
the approximate locations of the singularities of the functions involved. In 
particular, we have obtained our results without having to investigate the 
nature (or type) of the singularities that occur (e.g., we have not had to 
determine whether we have an algebraic or logarithmic singularity or a singu- 
larity of some other kind). In this sense, the methods we have used are 
quite powerful, especially in that all we really needed to know to make them 
work was a knowledge of the approximate location of the singularities, and 
not their exact location. It is interesting to point out, however, that the 
location of the singularities for this problem, being off the real axis, is 
in sharp contrast to the location of the singularities for "most" other 
physical problems which have been analyzed by similar methods (see e.g., Van 
Dyke [6]) . 

As far as the nature of the singularities is concerned, however, we have 
applied the method of analysis given by Hunter and Guerrieri [3] to the fre- 
quency series we developed in I and to the maximum amplitude series u(0;e), 

discussed in sections 3 and 4, which were both essentially series in powers 
2 

of e . The result of this analysis indicated that the singularities in 
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these cases are algebraic singularities with an exponent of +%, i.e., square 
root singularities. However, the method of [3] does not work well when more 
than one pair of complex conjugate singularities are about at equal distance 
from the origin, as in the case of the series (2.6) for x > 0, and hence we 
have no definite knowledge of the nature of the singularities for u when x > 0. 
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Figure Captions 

Caption 

Maximum amplitudes A(e) of the limit cycle calculated 
from (3.1) (dotted line) and the [48/48] Pade approximant 
(solid line) for 0 < e < 9 . The heavy dots are the 
values of amplitude computed using purely numerical 
methods by Zonnerfeld [7]. 

The zeros (dots) and poles (x's) in the upper half 

complex e-plane of the [48/48 ] Pade approximant to 

the amplitude series A(e) in (3.1). The distance of 

the singularities from the origin is denoted by R q 

and the angle of the singularity, measured from the 

positive real axis, is denoted by 8 . Here R =1.85 

o o 

and 8 = 0.8970, 

o 

Amplitudes A(e) computed using the transformed series 

(4.2) (dotted line), the series (4.4) (dashed line), 

and the asymptotic formula (4.3) (light dashed line), all 

expressed as functions of 6 given by (4.1), The 

numerical results of Zonnerfeld [7] are indicated by the 
heavy dashed line. 

Phase plane plots for the limit cycle u(x,t) computed 
using the first 100 terms of the series (2.6) with 
e = 0 (a), e = 1 (b) , e = 1.5 (c) , and e = 1.6 (d) . 

The location of the zeros (dots) and poles (x's) in 
the upper half complex e-plane of the [24/24] Fade 
approximant to the series (2.6) for x = tt/ 50 (a), 
x = 13TT/50 (b) , x = 18 tt/ 50 (c) , and x = 24 tt/ 50 (d) . 
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Captlon 

Schematic diagram of the upper half complex e-plane 
indicating the location of the fixed singularities 
(*) and the path taken by the moving singularity as 
x varies between 0 and it. The heavy dotted line in- 
dicates a portion of a circle about the origin, of 
radius R q , which corresponds to the radius of con- 
vergence of the series (2.6) at x = 0. 

Plots of the distance R = R(x) of the distance of the 
moving singularity from the origin (a) and the angle 
3 = 6(x) the singularity makes with the positive real 
axis (b) as functions of x for 0 < x < tt. 

Phase plane plots of the limit cycle of the Van der 
Pol oscillator using the transformed series (6.2), 
with w replaced by w given in (6.3), for e = 1.65 


(a) and e = 3.5 (b) . 
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Table II 

Comparison of the amplitudes computed from equations (4.2) and ( 4 . 4 ) 
with the numerical amplitudes reported by Zonnerfeld 7 . The percentage 

differences are enclosed in parentheses. 
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Table III 

Comparison of the amplitudes computed from equations (4.2) and (4.4) 
with the asymptotic formula (4.3) . The percentage differences are 
enclosed in parentheses. 
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